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A surface integral formulation is developed for the T matrix of a homogenous and 
isotropic particle of arbitrary shape, which employs scalar basis functions represented by 
the translation matrix elements of the vector spherical wave functions. The formulation 
begins with the volume integral equation for scattering by the particle, which is 
transformed so that the vector and dyadic components in the equation are replaced with 
associated dipole and multipole level scalar harmonic wave functions. The approach leads 
to a volume integral formulation for the T matrix, which can be extended, by the use of 
Green's identities, to the surface integral formulation. The result is shown to be equivalent 
to the traditional surface integral formulas based on the VSWF basis. 

© 2013 Elsevier Ltd. All rights reserved. 


1. Introduction and motivation 

The purpose of this note is to revisit the derivation of 
the transition (a.k.a., T) matrix for a homogeneous particle 
of arbitrary shape. The seminal work on this subject was 
performed by Waterman [1], the result of which was 
computational scheme that has various technical names, 
e.g., extended boundary condition method, null field 
method, yet which is commonly referred to, among the 
scattering community, as “the T matrix method”. Water- 
man's derivation begins with the vector Huygen’s princi- 
ple, which states that the exciting electric field inside the 
particle, and the scattered field outside the particle, can be 
related to the distribution of the tangential components of 
electric and magnetic fields on the surface of the particle. 
A representation of the surface fields in a vector spherical 
wave function (VSWF) basis ultimately leads to a formula 
for the T matrix which involves integrals, over the particle 
surface, of vector products of the VSWFs. Improvements in 
the computational scheme have been developed over the 
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years, to address highly aspherical particle shapes, chiral 
media, etc. [2,3]. 

An alternative method for T matrix calculations is via 
volume integral methods, such as the discrete dipole 
approximation (DDA) [4]. Recently, Litvinov demonstrated 
that, as would be expected, Waterman's surface integral T 
matrix formulation can be derived from the volume 
integral equation (VIE) formulation of Maxwell’s time 
harmonic equations [5], 

This paper will follow the same basic path as that cut 
by Litvinov, in that the analysis will begin with the VIE and 
will end with a surface integral T matrix formula. The main 
difference, however, is how the associated formulas are 
represented. In particular, the VIE will be transformed, at 
the outset, into a scalar form that employs the VSWF 
translation elements as basis functions. Subsequent opera- 
tions and transformations on the VIE, leading to the T 
matrix formulas, can then be done almost entirely in a 
scalar representation. 

2. Formulation 

Consider a particle of arbitrary shape, as illustrated in 
Fig. 1. The particle is taken to be homogeneous and 
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Fig. 1. System configuration. 


isotropic in composition, and characterized by a complex 
relative refractive index m. The medium in which the 
particle is immersed is taken to be non-absorbing. The 
interior and exterior regions of the particle are denoted as 
V int and V ext , respectively. Coupled with these volumes are 
two additional regions to be used in the formulation, 
which are defined by the circumscribing (radius r c ) and 
inscribing (radius r ( ) spheres centered about a fixed origin 
r 0 of the particle. The region V exti c is that external to the 
circumscribing sphere, and V intJ is internal to the inscrib- 
ing sphere. 

The electromagnetic scattering problem is typically 
described in a dynamical sense, i.e., a particle is excited 
by an incident field, which produces a scattered wave that 
propagates away from the particle. In a more mathemati- 
cally consistent description - which recognizes that the 
time variable is removed in the time harmonic formulation 
- the exciting field is the field which exists in the overall 
system when the particle is absent, and the problem is to 
predict the new field which results when the particle is 
present. The scattered field, in this context, is the differ- 
ence between the external fields with and without the 
particle present [6,7]. 

Since the exciting field contains no singularities in the 
vicinity of the particle, it can be described mathematically 
as an expansion in regular VSWFs centered about the 
particle origin. The scattered field, on the other hand, is 
described by an expansion in outgoing VSWFs about the 
particle origin, in order to automatically satisfy the far- 
field radiation condition [6]. The particle T matrix provides 
the relationship between the expansion coefficients for the 
exciting and scattered fields, so that 

oo n 2 

Esca(r) = 2 2 S ^mnpN ( m ^p(/<o(r— 1 * 0 )), reV eX f c (1) 

n = im = -np = i 

oo l 2 

a mnp = X X X T mn p klqf klq ( 2 ) 

1 = 1 k = -1 q = 1 


corresponding wavelength, and N' 1 ’ and N (3) denote the 
regular and outgoing VSWFs, respectively. These functions, 
and their associated properties, will be described in the 
subsequent section; for now the analysis will turn to 
the VIE. 

In the VIE formulation, the electric field at some point r 
is governed by 

E(r) = E e *c(r) + ko(m 2 -l) / £( /c 0 (r-r')) • E(r') dV (4) 

Jv, 


in which E exc (r) is the exciting field and g(k 0 (r-r')) is the 
free space dyadic Green's function, defined by 


g(k 0 (r-r')) 



exp(ik 0 |r-r'|) 
4zr|r— r' | 


(5) 


where ® denotes the dyadic product of two vectors, I is 
the identity dyadic, and i = \T- T. The dyadic Green's 
function has a singularity at r = r' which makes the actual 
solution of the VIE for the internal field nontrivial. This 
aspect of the VIE formalism is discussed by Van Bladel [8] 
and Lakhtakia and Mulholland [9] who show that when 
reV int the integrals in Eq. (4) should be taken as principal 
values around the point r'. The way we treat this implicit 
singularity will be clarified later. 

An operator-based solution to Eq. (4) can be con- 
structed by defining a dyadic transition operator T(r,r') 
so that, in regions within the particle [10,11], 

Eint( r ) = 7 i — \ / T(r,r') • E exc (r') d 3 r', reV int (6) 

ko(m 2 -l)7v int 


Replacing this into Eq. (4), and constraining r to lie within 
V int , gives 


L 


/<o(m 2 -l) 


nr, 


,r')-i5(r-r')I-^ Q(k 0 (r-r")) ■ T( r",r') d 3 r" 


•Eexc(r') d 3 r' = 0, reV int 


(7) 


This relation will identically hold providing T(r,r') satisfies 
r(r,r') = kg(m 2 -l)^(r-r')I + jT g(k 0 (r-r")) ■ T(r",r') dVj, 
r,r'eV int (8) 


and a general volume integral relation for the field at all 
points r can now appear as 


E(r) = E eX c(r)+ f g(k 0 (r-r")) • f T( r",r') 

J V int JVi nt 

•E ra c(r') dV d 3 r" 


(9) 


By using Eq. (6) in Eq. (4), the relationship between the 
scattered and exciting fields becomes 


Esca(r) — E(r)-E exc (r) 

= J S(k 0 (r-r')) ■ ^ T(r',r") ■ E exc (r") d 3 r" d 3 r', reV ext 

( 10 ) 


oo n 2 

Eexc(r) = X X X fmnpKnp( k o(r-r 0 )) (3) 

n=lm = -np=l 

where r is the position vector of the observation points, 
k 0 = 2 n/X is the wavenumber in the host medium, X is the 


Eq. (10) will provide the basic starting point to identifying 
the particle T matrix. Indeed, with some imagination, 
Eq. (10) can be seen as analogous to the VSWF T matrix 
relationships: the former involves a double integration 
over the particle volume, whereas the latter involves 
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a double summation (i.e., the row and the column order/ 
degree/mode indices) over the particle T matrix. Deriva- 
tion of the T matrix from the VIE transition dyadic will 
therefore involve the integral transformation of a spatial 
distribution into a VSWF expansion. Before proceeding 
down this path, however, it will be useful to introduce the 
needed features of the VSWFs. 


2.3. Vector spherical wave function background 


The VSWFs, of type t = 1 (regular) and t = 3 (outgoing), 
and order n, degree m, and mode p = 1 (TM) and p = 2 (TE), 
are defined here by 

/ 2 \ 1/2 

N ™ 2 ( kr > = (j^TPT)J V X (VmUkr)) HI) 


n1w=1vxn:;>) d2) 

where k is a positive real number and denotes the scalar 
wave function 


<n _ / Jn( kr )Y mn ( cos 9,<P), t = 1 
Vmn < r > - | h n (kr)Y m „( cos e,<p), t = 3 


(13) 


with j n and h n =j n + iy n denoting the spherical Bessel and 
Hankel functions and Y mn denoting the spherical harmonic 


y m „(cos 0,ip) = 


/2n + l (n-m)!\ 
V 4rr (n+m)\J 


1/2 

P™( cos P)e' m ^ 


(14) 


An essential component in transforming the VIE will be 
the ability to translate the basis functions from one coordi- 
nate origin to another. Translations of the VSWFs are 
performed by application of the addition theorem for VSWFs, 
which, for the problem at hand, will appear as [12,13] 


K'npOd r~r')) = 


oo 


Z 
/= 1 


/ 2 

2 2 Jklq mnp 

k = -lq = 1 


(-fcnhQfcr) 


(15) 


Nmnp(k(r-r')) 


= Z Z Z Jklq m „p(-kr')N®(kr), 

! = lfc = -lq = l 
X|r'| < in 


(16) 


Hmnp klq(k(T I* )) 

oo n' 2 

= 2 S 2 Jmnp m'n'p'(H(X~ ^o))^m'n'p' klq(H(V 0“^* ))» 

n' = lm' = -n'p' = 1 

x|r-r 0 | < |r 0 — r'| (19) 


Hmnp klq(H(V ** )) 

oo n' 2 

= Z Z Z fl/mip m „-p (k(r-r (l ))J m . n , p , k i q (k(r n -r )), 

n' - 1 m' = -n' p' = 1 

x|r— r 0 | > |r 0 — r'| (20) 

In the formulation developed here, the translation 
elements will take on a somewhat more elevated role, 
that being the set of basis functions for representing the 
fields in the VIE. This will have two distinct advantages in 
the formulation, being (1 ) the vector problem is reduced to 
a scalar problem, and (2) the basis functions become 
“transparently” translatable, as the components of the 
translation operation become the basis functions them- 
selves. An equivalence between the VSWFs and the trans- 
lation matrix elements can be obtained by letting |r| ->0 in 
Eqs. (15) and (17), for which the only surviving terms are 
those for the electric dipole (l=q= 1); this results in 

l'Cp(kr) = 2 i p iJku mnp(kr) (21 ) 


Nmnn(kr) = 


1 1 

n Z P k^kll mnp( k r) 


( 22 ) 


in which P fe is a unit magnitude Cartesian vector given by 


P_1 V2 


/1\ 


\oJ 


, Po = 


/ 0\ 

0 

W 



(23) 


The relation between the VSWF of different mode, i.e. 


u<i) 


— D.( kr >= l V x N mnnO<r) 


(24) 


oo / 2 

NmnpW-n) = z Z Z H«, m np(-kr')N«>(kr), 

l=lJc=— lq = l 

x|r'| > |r| (17) 

in which J and H are the regular and outgoing VSWF 
translation matrices; the elements of these identically satisfy 
the scalar Helmholtz equation and involve expansions of the 
type 1 and 3 scalar wave functions. The summation over 
order / in (15)— (17) formally appears as an infinite sum, yet it 
is understood that the series will converge to an arbitrarily 
small error in a finite number of orders, providing the 
convergence radii criteria are met. 

Eqs. (15)— (17) can be used to infer the following transla- 
tion and factorization properties of the translation matrices: 

Jmnp klq( k (^~r )) 

oo n' 2 

— Z Z Z Jmnp ^o))Jm'n'p’ klq( k (^0 ~ T )) 

n' = 1 m’ = -n' p‘ = 1 

( 18 ) 


will also lead to 

1 1 1 

Z Pm'/m'll mrn3-p)( kr l = ~T Z P m' 

m' = -1 K m' = -1 

x Z/m'll mnp( k T) (25) 

and the same relations would hold for the outgoing 
functions. The salient point here is that Eqs. (21) and 
(22) provide a Cartesian representation of the VSWFs. 

It can also be shown that the dyadic Green's function is 
related to the outgoing, dipole-level VSWFs by 

G(ka(r-r')) = 'V £ j 1 (k 0 (r-r' ) )<g>P* (26) 

Is! in m = -1 

where the asterisk denotes complex conjugate, and by 
combining the above with Eq. (22), the dyadic Green's 
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function becomes 

S(k 0 (r-r')) = ^ 2 Z U ku mll (/c 0 (r-r'))P k ®P* (27) 

D *m = -n=-1 


2.2. VSWF transformation of the VIE 

The VIE can be transformed so that its vector dimen- 
sions (i.e., the Cartesian components) are replaced with 
the -1, 0, and 1 azimuthal degrees for the dipole. This 
transformation of the transition dyadic produces a so- 
called two-point transition matrix T (2) (r,r'), in which 


TS(r.r') = ^P* • T(r.r') • P fc , m,k= -1,0,1 

(28) 

and since 


Pm ■ Pfc — fim.k 

(29) 

then 


T(r,r') = 2 Z TS(r,r')P m ®Pj 

K 0 m = -1 k = -1 

(30) 

and the internal and scattered fields, appearing in Eqs. (6) 
and (10), will now be given by 

Eint(r) = ’ 2 Pm/ 2 Or.r")Pj 
•E exc (r") d 3 r", rEV int 

(31) 

Esca(r) = X Pm/ 2 ^mll m'll(^o(J* — J*)) 

m = -1 J Vi nt m' = -1 


xf/ z C(r'.r")P: • E exc( r") d 3 r"' 
\Jv int k = - 1 ) 

r<=V ext 

) d 3 r, 

(32) 

in which 


a = ^(m 2 -l) 
on 

(33) 


Using the definition in Eq. (30) along with the VIE of 
Eq. (8) results in a VIE for T <2) : 

-r^( r , r ') = 5(r- r ')5 m _ k 

a mK 

+ [ 2 H mll m 'ii(fc 0 (r-r"))T® fc (r",r') d 3 r", 

■J V int m' = - 1 

r,r'EV int (34) 

It should be emphasized that Eq. (34) is fundamentally 
equivalent to Eq. (8); the only distinction is that the former 
works in dipole degree space, whereas the latter works in 
vector component space. The wave function representa- 
tion, however, does make it possible to define additional T 
matrix operators, that arise once Eq. (10) is formally 
integrated over one or both of the volume domains. These 
operators will involve wave function expansions beyond 
the dipole level, and as such they will have no equivalence 
to operators in the vector component model. 

Performing such integrations, however, requires that 
certain restrictions be made on the nature of the exciting 
field and the location of the evaluation point r in Eq. (10). 


Specifically, the exciting field is taken to arise from sources 
located in V extiC , i.e., points outside the circumscribing 
radius centered about the particle origin r 0 . A sufficiently 
general model representation of the exciting field, for this 
restriction, is to take it as that radiated from a single 
dipole, of some specified orientation, located in V exCC ', an 
arbitrary exciting field - including a plane wave - could be 
constructed from a superposition of dipole sources. 
Assume that the exciting field originates from a point 
dipole source located at position T e sV extiC . The orientation 
of the dipole is specified by unit vector u d relative to the 
particle coordinate system, and the field radiated by the 
source can be described by 

P*.E exc (r) = j; H mll M1 (ko(r-r e ))s fc , r<=V extiC (35) 

where E 0 is a characteristic electric field amplitude, the 
presence of which renders the s k coefficients dimension- 
less, and 

Sk = Pk • u rf (36) 

Since the dipole source is located outside the circumscrib- 
ing sphere, Eq. (19) can be used to represent the exciting 
field at points internal to the particle as a regular VSWF 
expansion, given by 

P* • E exc (r) = -|L Z/mii „(ko(r-r 0 )) 2 M i(ko(r 0 -!•«■))<>* 

A/ fc=-l 

= ^J=IWM r ~r 0 ))/°. reV int (37) 

V J7Z H 

In the above and what follows, Greek subscripts will be 
shorthand for the triplet of degree/order/mode, i.e., 
p = (mnp); this convention will greatly simplify the result- 
ing presentation. 

The general form of Eq. (37) will hold for any dipole 
source located in V ext c , in that alterations of the exciting 
field would only affect the expansion coefficients f°. More 
importantly, the form enables the formal integration over 
r" in Eqs. (31) and (32), the process of which defines a one- 
point transition matrix Tj6(r',r 0 ), so that 

Pm-Eint(r) = £ " ir^lr.ro)/! 1 , reV int (38) 

av in v 

Pm-Esca(r)= ~^= 2 l/ H m ii m .n(kb(r-r'))7^? e 

V in m' = — 1 v JVj nt 

(r',r 0 ) d 3 r'/°, reV ext (39) 

where T (1) is defined as 

r m>' r o)= Z / T m k ( r ’ r ')Jku Jk 0 (r'-r 0 )) d 3 r' (40) 

k = -1 Jv int 

A VIE for T (1) is obtained by multiplying Eq. (34) into the 
regular translation matrix y /(11 ^(kolr'-ro)) and integrating 
over r', to yield 

l T ml( r,r 0 ) =J m n „(ko(r-r 0 )) 

+ f 2 H m n m-nlMr-r'pr^^r'ro) d 3 r', 

JV inc m- = -1 

reVint 


(41) 
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The same restriction that was placed on the exciting 
field can now be applied to the scattered field. Specifically, 
the observation point r of the scattered field is constrained 
to be located solely in V e xt,c- Under this restraint the 
outgoing translation matrix H(k 0 ( r— r')) can be factored 
per Eq. (20). And as before, this allows for the integration 
over interior points r' to be separated from the location of 
the evaluation point r. The end result is the T matrix 
relationship 

P* ,E sca (r) = £H m „ ^(/co(r-r 0 ))2T p „(r, 0 )/°, reV ext , c 

V n v 

(42) 

in which 

r„„( r 0 )= 2 [ y„ mll (fc 0 (r 0 -r'))T^(r',r 0 )d 3 r' (43) 

m = -1 JVj nt 

T„„(r 0 ) = 2 Z / J H m n(ko(r 0 -r')) 

m = -1 k = -1 -'Vint 

*f v ^i(r',r"l/kn „(fe 0 (r"-r 0 )) d 3 r" d 3 r' (44) 

The integral relations for the matrix operators T^fr.ro) 
and T,, „(r 0 ), given in Eqs. (41) and (43), are equivalent to 
those derived by Litvinov [5], They are also analogous to 
the formulas used to construct the T matrix for a cluster of 
spheres, via the superposition T matrix method [12], And 
along the same lines as the superposition method, Eq. (41 ) 
provides an energy conservation statement for the T 
matrix. This is obtained by multiplying the equation 
through by T^*\ r,r 0 ), summing over m and integrating 
over r, and employing the properties of the translation 
matrices for real k 0 : 

J^fkoir-To)) =J* l (k 0 (T 0 -r)) (45) 

H m u m(ko(r-r'))+H*„ mll (k 0 (r'-r)) = 2J mn ku (k 0 ( r-r')) 

= 2 U* mu (ko(r 0 -r ))J H M i(fco(r 0 -r')) (46) 

b 

This results in 

-Ref-) Z [ l^(r,r 0 )| 2 d 3 r+2|T^(ro)| 2 = -ReT TO (ro) 
\ a / m = - 1 Jv int e 

(47) 

When summed over v and multiplied by 7t/kg, the two 
terms on the left correspond to the random orientation 
absorption and scattering cross sections of the particle; 
note that Re(l /a) will be zero for non-absorbing material. 
These equate to the random orientation extinction cross 
section, as given by the right term. 

3. Calculation of the T matrix 

3.1. Surface integral relations 

Eq. (41) provides the starting point to calculate the T 
matrix for a set particle geometry and composition. The 
most direct implementation, in this regard, is to adopt a 
discretized volume integral formulation, for which the 
integration over volume is replaced with a sum over 
discrete control volume elements, and the quantity a is 


replaced by an effective polarizability of the element 
[14,15]. The result is basically the DDA for homogeneous 
particles, applied to Tj^(r;,ro) at discrete cell positions r,- 
and for a set degree/order/mode u component. Upon 
solution for a given v, the contribution to T^(r 0 ) would 
be obtained by the discretized volume integration in 
Eq. (43). The practical implementation of this approach 
is discussed in [4]. 

On the other hand, a formal integration of the volume 
integral in Eq. (41) will lead to a T matrix calculation 
scheme based on a surface integral. The surface integral is 
derived from Green's second identity and by recognizing 
that the translation matrix H and the one-point matrix T (1) 
satisfy the scalar Helmholtz equation (SHE) with wave- 
numbers k 0 and mfc 0 (the latter comes from the association 
of T (1) (r,r 0 ) with the internal field via Eq. (38)). This 
procedure will split the volume integral into two surface 
integrals; one over the particle boundary, and another over 
a vanishingly small spherical surface centered about the 
singular point at r' ->r, plus a third part which accounts for 
the singular point contribution [8,9]. To evaluate the 
second surface integral, we use the fact that T (1, (r',r 0 ) 
must contain no singularities for r'eV int , and that it must 
be locally expandable in regular scalar spherical wave 
functions - or, equivalently, in regular translation matrix 
elements. Consequently 

lim 7 m!( r '. r o) = Z JmU mOl(kl( r '- r ))^nr iX r ,r 0 ) (48) 

r ^ r m' = - 1 

and the VIE of Eq. (41 ) becomes 
l T ml< r ’ r a)=Jmu „(fco(r-r 0 )) 

+ ^T-r Z [Q[H m iini(ko(r-r')),T^ ) (r',r 0 )] r ., B 

m k = - 1 

-Q[H m „ ku(-k 0 r')J k u nui(kir')] r -,s„^(r,r 0 )] 

r<=V im (49) 

K 

where B and S 0 denote the particle surface and a spherical 
surface of vanishingly small radius, and the surface inte- 
gral operator Q is defined as 

Q[F(r),C(r)] rS =X [n- (VF(r)G(r)-F(r)VG(r)) d 2 r (50) 

k 2 Js 

with n and S denoting the outward surface normal and the 
closed surface of integration, respectively. 

For a spherical surface - such as that enclosing the 
singular point in Eq. (49) - the surface integral can be 
analytically evaluated to yield (in the general case) 

1 27T\ t 9 

Z Q[H,, mil (ka(— r)), J mU „(/q,(r))] r So = 2 3 (2k a +k b )S llIJ 
m = -l k 0 k a 

(51) 

Algebra will then reveal that the sum of the second and 
third terms on the right hand side of Eq. (49) will 
identically cancel the left hand side for all reV,- nt . 
This results in 

Jmu „(ko(r-r 0 )) 


118 


D.W. Mackowski et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 123 (2013) 113-121 


= , 2 Q[H m , , k , , (fco(r-r')),7^ ) (r',r 0 )] r , B , 

re Vint (52) 

Eq. (52) is equivalent to the vector Huygen’s condition for 
points inside the particle, i.e., the null-field condition [1,3], 
This can be seen by multiplying through by Pmf v and 
summing over m and v , for which the left hand side 
becomes the exciting field at interior point r and the right 
involves a distribution of tangential interior surface fields. 

The volume integral relation for the T matrix, Eq. (43), 
can be put in a surface integral form by the same 
procedure, leading to 

f^( r o) = -J—r 2 Q[/^n(feo(ro-r')),T^(r',r 0 )] r . iB (53) 

m -* k = -1 


3.2. Representation of the interior fields 

The surface integral relations in Eqs. (52) and (53) are an 
alternative form to the vector relations initially derived by 
Waterman [1]. Getting to this point has basically involved a 
contraction of the dependent variables in the VIE. Specifi- 
cally, the sought quantities in the surface integral relation of 
Eq. (52) are T^*(r,r 0 ) and n ■ VT^’(r,r 0 ) at surface points 
reB, and these quantities are sufficient to calculate the T 
matrix via Eq. (53). Indeed, it can be easily shown, via a 
volume-to-surface integral conversion of Eq. (39), that field 
at all external points can be obtained from a surface integral 
relation analogous to Eq. (53). 

The interior field (by which we mean T5^(r,r 0 )) can also 
be obtained from a surface integral relation, although one 
that is not derived from the VIE. The principles applied to 
derive this relation are, again, based on the facts that 
T^(r,r 0 ) obeys the SHE, and that it contains no singula- 
rities for reV int . Because of this, Green's second identity 
will state that 

2 Q[Hmii m'ii(fei(r-r')),T ( n J>„(r',ro)] r - iS(I . ) = C (54) 

m' = -1 

must evaluate to the same value C for all closed surfaces 
S(r) enclosing the point r. This surface can be shrunk to a 
vanishingly small sphere centered on r - for which 
Eqs. (48) and (51) can be applied to analytically evaluate 
the integral - and it can also be extended to the interior 
surface of the particle. This results in 

O r , r o) = 2 Q[W m i i mu i (fci (r-r'» , r^(r ' , r 0 )] r . B , 

fc 0 m m' = -1 

reVfnt (55) 

The next step in the problem has to do with the 
representation of the field in an analytical form; such a 
form will be needed to close the problem and obtain 
practical matrix equations for T. An application of the 
translation in Eq. (19) to Eq. (55) would result in the 
regular expansion 

f ( m'(r,r 0 ) = XJ mu „(ki(r-r 0 ))W„„(ro) (56) 


with the coefficient matrix W given by 
W„„(ro)=^ X QlHp m'ii(ki (r 0 -r')),T^(r',r 0 )] r , iB 

m' = -1 

(57) 

yet, per the restriction on Eq. (19), this expansion would 
appear to be valid solely for points r within the inscribing 
sphere centered about r 0 . How valid is it, then, to extend 
the expansion in Eq. (56) to the surface? 

The short answer to this question is yes: Eq. (56) can 
represent, in a least-square-error sense, the surface fields. 
However, providing a convincing yet easily accessible 
explanation as to why this is the case continues to be, in 
our opinion, one of the most challenging and frustrating 
aspects of the entire effort. Indeed, Waterman's 1979 
reformulation of his original 1971 derivation was largely 
focussed on the issue of interior field representation and 
the validity of the so-called internal Rayleigh hypothesis 
(that being the use of Eq. (56) throughout the interior 
volume) [1,16]. We will attempt to offer a simple rationale 
for the use of Eq. (56), and will refer the reader to Ref. [3] 
for a comprehensive mathematical proof of the complete- 
ness of the regular wave functions. 

The most salient point, to us, is that an expansion of the 
field in regular wave functions - such as Eq. (56) - will 
automatically satisfy Eq. (55) on a term-by-term basis. 
That is, application of the surface integral principle behind 
Eq. (55) to the regular wave function will yield 

Jmi i „(ki(r-r 0 )) 

1 

= 2 Q[Hmii m'ii(ki(r-r'))J m ii p(ki(r'-r 0 ))] r , iB , reV int 

m' = - 1 

(58) 

This result, it should be emphasized, has nothing to do 
with the translation theorem. Now multiply this relation 
through by the coefficient matrix W^(r 0 ), defined in 
Eq. (57), restrict r to points inside the inscribing sphere, 
sum over order/degree/mode p up to some order trunca- 
tion limit L, and subtract the result from Eq. (55). Denoting 
the approximation error at point r, et(r), as the difference 
between the actual field value T^’(r,r 0 ) and the truncated 
regular wave function expansion evaluated at r, this 
exercise results in 

|%(r)= 2 Q[H m ii m ii(ki(r-r')),e L (r')] r . B , reV inti , 

k 0 m m' = - 1 

(59) 

Since points within the inscribing sphere meet the con- 
straints on the translation theorem for casting Eq. (55) into 
Eq. (56), the error on the left hand side can be made 
arbitrarily small by increasing L. The limit of 6i(r)->0 can 
only be satisfied by two surface distributions in Eq. (59): 
et(r')->0 or er(r')->H m ll mll (ki(r'-r)). Yet this second case 
cannot be present as the internal field contains no singu- 
larities. The implication is that Eq. (56) can approximate, to 
an arbitrary precision depending on the truncation limit, 
the internal field at all points— including the surface. 
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3.3. The matrix equations 

When Eq. (56) is substituted into Eq. (52), one obtains 
Jmu ,(Mr-r 0 )) 

= - -J-f S I QIHmii kuO<o(r-r')),J ku ,(/q(r'-r 0 ))] r , B 

m — 1 fc = -i 

»(r 0 ), rGV in t ( 60 ) 

It is assumed that the expansion for the surface fields is 
truncated at some set order I, and this, in general, would 
result in 2L(L+2) unknown coefficients in W for a set 
column order/degree/mode state u. The final task of the 
problem is generating a system of equations for the 
coefficients. 

The approach followed by Waterman was to restrict r to 
points inside the inscribing sphere, for which the transla- 
tion theorem in Eq. (19) can be applied to cast the right 
hand side as expansion of regular wave functions centered 
about r 0 . Since the left hand side is also a regular wave 
function, this will imply that W must satisfy 

- V = IQ'I' 1, (r 0 )W, „(r 0 ) (61 ) 

v' 

with 

Q B 3,1) = I Q[H b tii(ko(r 0 -r'))J M i „(mko(r'-r 0 ))] r , iB (62) 

k = - 1 

The constant factor fco(m 2 -l) has been absorbed into W in 
going from Eqs. (60) to (61); this quantity would cancel 
out in the end. The T matrix is obtained from Eq. (53), 
which would now appear as 

r„„( r 0 ) = „(r 0 ) (63) 

v' 

where Q 1,1 is obtained from Eq. (62) with H replaced by J. 

It is instructive to examine how the derived formulas 
remain consistent to the energy conservation requirement 
of Eq. (47). By using the translation properties of regular 
wave functions in Eq. (18), it can be shown that 

SlV(r 0 )l 2 = I,(Q!^\r 0 )WMr 0 ))(Q})'}\r 0 )W u . „(r 0 ))* 

= QtVrnu „-(ki( r-r 0 ))VW Jr 0 ))*, 

Oil mu im(ko(r-r'))Jkii „.(ki(r'-r 0 ))W„. „(r 0 )] r , B ] r , B (64) 

in the above and the following equation, summation over 
all subscripts not appearing on the left hand side is 
implied (i.e., summation over everything except v). The 
double surface integral can be converted directly to a 
double volume integral since the integrand is free of 
singular points. And since the volume integral over r 
would be exactly the conjugate of the volume integral 
over r', the regular wave function J(k 0 ( r-r')) can be directly 
switched with H(k 0 (r-r')); this is basically going back- 
wards from the derivation of Eq. (47). The volume integral 
can now be converted back into a surface integral, yet this 
will now require accounting for the singularity in H per the 
procedure in Eq. (49). This results in 

D|T^(r 0 )| 2 = Q[(/ mll „*(ki(r-ro))W„. „(r 0 ))*, 

QIHmii kii(ko(r-r'))Jkn „-(ki(r'-r, 0 ))W„. ^(r 0 )] r , iBo ] r B . 


+ Re 


T 

a 


QlUmn A k i(r-r 0 ))W„. „(r 0 ))*, 


C/mii „'(ki(r-ro))W„- „(r 0 ))*] r , B 


(65) 


Both of the surfaces B 0 and B t in the above relation are 
formally on the particle surface; the only distinction is that 
Bj is enclosed by B 0 . That is, B 0 contains the singularity in 
H, yet B , does not. Eqs. (55) and (63) can be now applied to 
show that the double surface integral reduces to 
-Re T„„( r 0 ), and the last term is simple the surface integral 
form of the absorption contribution in Eq. (47). 

As would be expected, the surface integral formulas for 
T are shown to give a result that - for a sufficiently large 
truncation order - is consistent with energy conservation. 
The caveat with respect to the truncation order is impor- 
tant here, as the derivation relied on formally taking the 
number of row orders in Eq. (64) to infinity in order to 
reduce the product of two regular translation matrices to a 
single matrix per Eq. (18). Waterman recognized that his 
method did not automatically satisfy energy conservation 
for a finite truncation order, and he devised a means of 
enforcing the unitarity property of T (for Re[fcj]) via an 
orthogonalization of the matrix equations [1], It is inter- 
esting to note that Eq. (65) suggests an alternative matrix 
formulation which - potentially - would automatically 
satisfy energy conservation at each order, that being 

“(0*= £ ZQ[C„„(k iC-ro», 

m,k = -1 //' 

Q[W m n fci i (fco( r — r '))i/ ki i ^(/q(r -ro))]r, Brj ]r.B, 

W„ „(r 0 ) (66) 


Per the discussion with regard to the internal field, the 
function H miuil (fc 0 (r-r')) could be represented as an 
expansion of regular and outgoing parts, and this would 
split the coupled double integral into a product of two 
integrals - and doing so would simply recover Eq. (61). 
We do not know if implementation of Eq. (66) is practical; 
certainly, one would have to devise some way of dealing 
with the near-singular point at r = r'. 


4. Discussion 

It occurred to the lead author, at the completion of the 
derivation leading to Eqs. (61) and (63), that a more direct 
route might have been to simply substitute the Cartesian- 
based formulas for the VSWFs, given in Eqs. (22) and (21), 
into Waterman's surface integral formulas for the Q. and 
RgQ matrices. This is done now to demonstrate the 
veracity of the derivation. 

The conventional formula for the T matrix can be 
written, in usual matrix notation, as [1,3] 

T = -Rg Q. Q _1 (67) 

where (assuming nonmagnetic media) 

= J n ■ ((V x NL 3, (/< 0 (r-r 0 ))) x N' l) (mk 0 (r-r 0 )) 

+ N( 3, (ko(r-r 0 )) x (V x N^'(mMr-r 0 )))) d 2 r (68) 

with denoting the triplet (- mrtp ) for / 1 = ( mnp ), i.e., 
switched signs on the azimuth degree. The formula for 
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Rg Q. is the same, with the outgoing VSWFs replaced by 
the regular VSWFs. 

Denote the vector part of the integrand in Eq. (68) that 
is dotted into the normal as K^. In terms of the Cartesian 
based VSWF formulas in Eqs. (21) and (22), this quantity 
would appear as 

= n I I ((V X „(r» X Pfc g k , „(r) 

- >n m' = -1 k' = -1 


+ Pm'/ m> (r) X (V X P k g k , „(r))) 

(69) 

where the scalar functions / and g are 


/ m >(r) = H m .„ ? (Zc 0 (r-r 0 )) 

(70) 

gk ■ t/( r ) =7^11 „(mk 0 (r-r 0 )) 

(71) 

By using the vector calculus identities 


V x (a/(D) = V/(r) x a 

(72) 

a x (b x c) = b(a ■ c)-c(a ■ b) 

(73) 


and switching the m' and k' indices to group terms, the 
integrand can be written as 

= — q— 2 2 ((V/ m . „(r)g k , Jr) 

~ >7Z m' = -1 k' = -1 V 

-fm' „(T)Vg k , „(r)(P m . ■ P*0 

-p m [Pk- • (V/ m , „(r )g k , „(r )-Vg m , „(!•)/„, „(r»]) (74) 

The P vectors have the property 

Pm'Pk'=( — f) />m',-k' (75) 

which, when applied to the first term in Eq. (74), gives a 
result that can be reduced, in a few steps that make use of 
the symmetry properties of the translation matrices, to the 
form of the integrand in Eq. (50). 

One would therefore expect that the second term in 
Eq. (74) must be equivalent to the first, or that it must be 
zero. The actual situation is somewhat more complicated: 
the second term is not zero, and it cannot be reduced to 
the form of the integrand in Eq. (61 ), yet its net contribu- 
tion to the T matrix will be zero. To show how such is the 
case, the term is expanded via 

2 2 Pm'[P*' ■ (V/ m> (r)g k , „(r)-Vg m , „(r)/ k> (r))] 

m' = -1 k' = -1 

= 22 Pm'[Pk'-V(f m> (r)g lk , y (r)-g m ,„(rl/ 7 Jt ,„(r))] 

m! = -1 k' = -1 

- 2 2 Pm'[Pk' • (V& „(r)/ m> (r)-V/ fc> (r)g m . ^(r))] 

m' = -1 k' = -1 

(76) 

The VSWFs have zero divergence, which implies that 
2 V-(P k ./ k> (r))= £ p r v/ t> (r) = 0 (77) 

k • = -1 k 1 = -1 

and likewise for the g function; this property eliminates 
the last term in Eq. (76). The second term in Eq. (76) is not 
zero, yet its contribution to the surface integral can be 
evaluated by using the divergence theorem: 

2 £ f n.(P m -[Pk'.V(/ m . ^(r)g k , „(r) 

m' = -1 k' = -1 JB 


-gm> AVfk' „(r))]) d 2 r 

= ££/ V-(P m .[P*. V(fm H (T)g k , v (t) 

m' = -1 k' = -1 •' V i n f 

-gm’ J T )fk’ „(*•))]) d 3 r 

= 22 / Pm' • V[P k . ■ V(f m , H (r)g k , „(r) 

m! = -1 k' = -1 J V int 

-gm'v(r)fk'„(r))]d 3 r (78) 

The order of the P ■ V operators, and the m' and t indices, 
can be switched in the second term on the last line, to 
show that the integrand in the volume integral is zero. 

An additional check on the formulation can be obtained 
by applying Eqs. (61) and (63) to a spherical particle. It has 
been verified by the authors that the T matrix for this case 
is diagonal and azimuth degree degenerate, and with 
elements equal to the Mie coefficients. 

In closing, it has been demonstrated that the surface 
integral formulas for the T matrix of a homogeneous, 
isotropic particle can be formulated entirely with the set 
of scalar basis functions provided by the VSWF translation 
matrix elements. No claim is made that the resulting 
formulas offer any numerical advantage over the tradi- 
tional, VSWF-based formulation. However, the new for- 
mulation could, conceivably, result in some coding 
simplifications. 
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Appendix 

Explicit formulas for the translation matrix elements 
that appear in this work are 


Jm 11 fc/l (^ <r ) — 


6n 




1/2 


(2/-1 )(2Z+ 1)(2Z +3)/ 

(f+l)(2f+3)((-k)!(f+k)! 

!(l-m)!(m+l)!(Z+fc-m-l)!(Z-k+m-l)! 


1/2 


x v4-m l-jO <r > 

,m / Z(2Z— l)(Z+fc— m + 1 )!(Z— Zc+m + 1)! \ 1/<2 
V (Z+l)(l-m)!(m+l)!(f-k)!(Z+k)! ) 

x^L^kr)] (79) 


Jm 11 ka0 <T ) 



= i(m(Z+l)-Zc) 

6 jr(l-k)\(k+l)\ 

l)(i-m)!(m+-l)!(l+k-m)!(l-k+m)! 


1/2 

¥4-m ,('«> 


( 80 ) 
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The formulas for the outgoing matrix H are the same, with the 
regular scalar wave harmonics replaced by the outgoing type. 

The gradient of a scalar wave function is given by the 
formulas 


+ '-^V mn (kT) 


= k 


1 


2n+l 


1/2 


(n-m-l)(n-m)\ 1/2 


2n-l 


V^m + l n 1 1 kv) 


(n+m + l)(n + m + 2)\ 1/2 
2n^3 ) V/m+ '"+' (kr) 


(81) 


d . d 

1 — 

Ox dy 

= -k 

+ 


Vmn( k r) 

1 X '/2 


2n + l 


(n + m-l)(n + iTi)\ 1/2 

2 n-i ) 

(n-m + l)(n-m+ 2)^ 1/2 
2h+3 ) ^m-ln + l(fa) 


V^m-1 n-l(fer) 


(82) 


dz 


VmnO < r ) = k 


1 


2n + l 


1/2 


(n-m)(n + m)\ 1/2 

2n-l J 


Y'mn-i(kr) 


(n-m+l)(n+in + l)y /2 
2n+3 


Vimn+i(fcr) 


(83) 
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